source("rperm_fun.R")

args <- commandArgs(TRUE)
j <- eval( parse(text=args[1]) )
output <- paste("did-", j, ".txt", sep="")

set.seed(48109)

q1 <- 6
q0 <- 6
t0 <- 10
t1 <- 10
rho <- .5
gamma <- .8
beta <- c(1,1,1)
theta <- 1
n.deltas <- 4
deltas <- seq(0, 3, length.out = n.deltas)

sims <- 10000

idx <- MakeIndex(q1,q0)
q <- q0 + q1
n <- t0 + t1
brn <- 500
nbrn <- brn + n 

#j <- 1 # 1, 2 , 4
sigv  <- c(rep(1, q - j), rep(20, j))
sigw  <- c(rep(1, q - j), rep(20, j))
sigx1 <- c(rep(1, q - j), rep(20, j))
sigx2 <- c(rep(1, q - j), rep(20, j))
sigx3 <- c(rep(1, q - j), rep(20, j))

id <- rep(1:q, rep(n, q))
treat <- rep(0, q*n)
treat[1:(q1*n)] <- 1
iv <- rep(c(rep(0, t0), rep(1, t1)), q)
ivtreat <- iv*treat

totres <- matrix(NA, n.deltas, 4)

for(h in 1:n.deltas) {

	res <- matrix(NA, sims, 4)
	delta <- deltas[h]

	for(i in 1:sims) {
		u <- matrix(NA, nbrn, q)
		u[1, ] <- 0
		for(t in 1:(nbrn-1)) {
			v <- rnorm(q, 0, sigv)
			u[t+1, ] <- rho*u[t, ] + v
		}
		u <- u[(brn+1):nbrn, ]
		u <- c(u)
	
		x1 <- matrix(NA, n, q)
		x2 <- matrix(NA, n, q)
		x3 <- matrix(NA, n, q)
		mat.ivtreat <- matrix(ivtreat, n, q)
		for(k in 1:q) { 
			x1[, k] <- gamma*mat.ivtreat[, k] + rnorm(n, 0, sigx1[k])
			x2[, k] <- rnorm(n, 0, sigx2[k])
			x3[, k] <- rnorm(n, 0, sigx3[k])
			}
		x1 <- c(x1)
		x2 <- c(x2)
		x3 <- c(x3)
		
		y <- 1 + theta*iv + delta*ivtreat + beta[1]*x1 + beta[2]*x2 + beta[3]*x3 + u
		
		d <- data.frame(id, y, iv, treat, ivtreat, x1, x2, x3)
		
		f0 <- lm(y ~ iv + ivtreat + x1 + x2 + x3, data = d)
		#summary(f0)
		tstat <- coef(summary(f0, cluster = "id"))[3, 3]
		restric.f <- lm(y ~ iv + x1 + x2 + x3, data = d)
		restric.resid <- residuals(restric.f)
		restric.fit <- fitted(restric.f)
		boot.t <- NA
		for (b in 1:199) {
			d$boot.y <- restric.fit + rep(sample(c(1, -1), q, replace = TRUE), rep(n,q)) * restric.resid
			boot.t[b] <- coefficients(summary(lm(boot.y ~ iv + ivtreat + x1 + x2 + x3, data = d), cluster = "id"))[3, 3]
		}
		
		that <- rep(NA, q)
		for(k in 1:q) {
			that[k] <- coef(lm(y ~ iv + x1 + x2 + x3, data = d[d$id == k, ] ))[2]
		}
		res[i, 1] <- as.numeric(tstat > qt(.95, q-1))
		res[i, 2] <- IMTest(that, .05, q1, q0)
		res[i, 3] <- OrdTest(that, 903, idx, q1, q0)
		res[i, 4] <- as.numeric(tstat > quantile(boot.t, .95, type = 1))
	}
	
	totres[h, ] <- colMeans(res)
	write.table(totres, output)
}
